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Abstract. We present hydrodynamical simulations, using a 2-D two component model (ambient medium and 
pulsar wind have different specific heat ratios), of bow shocks in a representative regime for pulsar wind driven 
bow-shock nebulae. We also investigate the behaviour of a passive toroidal magnetic field wound around the pulsar 
velocity direction. Moreover we estimate the opacity of the bow-shock to penetration of ISM neutral hydrogen: 
this quantity affects observable properties of the nebula, like its size, shape, velocity and surface brightness 
distribution. Finally we compare these numerical results with those from an analytical model. The development 
of more realistic models is needed in order to tune the criteria for searches of new such objects, as well as to 
interpret data on the known objects. 
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1. Introduction 

Pulsars moving supersonically with respect to the ambi- 
ent medium are expected to give rise to a bow shock. 
In the case of interaction with the interstellar medium 
(ISM) the Ha emission from the nebula may be detected. 
Although pulsars have high typical velocities, and there 
are more than one thousand known radio pulsars, only 
four bow shocks have been discovered so far: PSR 1957+20 
(Kulkarni & Hester |198SD , PS R 222 4-f65 (Cordes et al. 
199^ ), PSR J 0437+ 4715 (BeU |1995D , and PSR 0740-28 



(Jones et al. 2001); while it is not clear whether the 
cometary-like nebula associated with the isolated neu- 
tron star RX J185635~3754 (van Kerkwijk & Kulkarni 



2001) is a bow shock or results from photo-ionization. In 
this kind of object the relativistic particles in the pul- 
sar wind produce a synchrotron emission too weak to be 
detectable in radio wavelengths (Gaensler et al. 20001 ). 
These nebulae may be detected instead in optical Balmer 
lines (mostly Ha ) as a signature of a non-radiative shock 
moving through a partially neutral medium (Chevalier 
& Raymond 1980). The Balmer lines are due to the de- 



excitations of neutral H atoms, following collisional exci- 
tations or exciting charge-exchange processes. When re- 
combination times are long these are the dominant pro- 
cesses giving emission in the optical band, as it has been 
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widely studied on various supernova remnants (Smith et 
al [1991D . 

A puzzling point is that two out of the four known neb- 
ulae have a shape closely matching that of a classical bow 
shock, while the others show a more peculiar shape with 
a conical tail: the nebula associated with PSR 2224-1-65 
in fact shows a remarkably conical tip followed by a bub- 
ble (which justifies why it has been named the "Guitar 
Nebula"); (Cordes 1996 ). Various hypothesis have been 
put forward to justify its shape: a peculiar ISM density dis- 
tribution in the surroundings of the PSR 2224-1-65 (Cordes 
et al. 1993); effects like mass loading due to neutral atoms 
which might penetrate in the external layers of shocked 
material (Bucciantini & Bandiera 2001, hereafter Paper 
I) . The fourth nebula, which has recently been discov- 
ered (Jones et al. |2001 ) near PSR 0740-28, seems to rep- 
resent an intermediate case, with a "standard" head and 
a conical tail. 

A pulsar wind interacting with a homogeneous ISM 
is in principle a much easier problem to model than that 
with a surrounding supernova remnant, since in the latter 
case the ambient medium might be quite inhomogeneous. 
Thus our analysis will be restricted to the former case. The 
ISM is seen by the pulsar as a plane-parallel flow, likely 
with a constant density (at least on the typical length 
scale of the bow shock) , and lasting long enough to pro- 
duce a steady-state regime. However, modeling is made 
more complicated by the fact that the pulsar wind is rel- 
ativistic and magnetized. Moreover the ISM is typically 
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partly neutral and the neutral atoms may have coUisional 
mean free paths comparable to the scale length of the sys- 
tem, then invalidating a purely fluid treatment. However 
as demonstrated in a previous paper (Paper I), if we sup- 
pose that the H atoms can be ionized only via collisions 
with electrons and protons of the shocked plasma, for a 
large number of pulsars the presence of a neutral com- 
ponent in the ISM can be neglected as far as the fluid 
dynamics is concerned, or at least it can be taken into 
account as a small perturbative effect. 

The interaction of the relativistic magnetized pulsar 
wind with the ISM ionized component is not direct, be- 
cause the relativistic particles themselves have very small 
cross sections, but is mediated by the magnetic field ad- 
vected by the pulsar wind and compressed on the head 
of the nebula. The effective mean free path of particles 
is then their gyroradius, which is typically much smaller 
than the typical dimension of the nebula, scaling with the 
distance of the bow shock stagnation point from the pul- 



(1) 



where C, and are the pulsar luminosity and peculiar 
velocity, and the density of the dragged component 
of ambient medium (Paper I). For this reason we have 
assumed little mixing between the two media, which allows 
the use of a fluid approach to investigate the structure and 
dynamics of these objects. 

In Section 2 we present the numerical code and the 
parameters of the various simulations we have performed; 
Sect. 3 discusses the results of these simulations and com- 
pares them to what is expected fr om an analytic "two thin 
layer" model (Comeron & Kaper 1998 ); in Sect. 4 we eval- 
uate the penetration thickness of the external layer for a 
neutral hydrogen atom of the ISM. 

2. Simulation Setup and Initial Conditions 

To simulate a spherically symmetric wind from a moving 
star, interacting with the ambient medium, we have used 
the multi-D code CLAWP ACK (LeVeque |l994[ ): a second 
order Godunov (Godunov |l959 ) with piecewise linear re- 
construction and transverse flux corrections. The problem 
has an intrinsic cylindrical symmetry which allows us to 
use a two-dimensional grid with Strang Splitting (Strang 



196S) for the geometrical terms. The boundary conditions 



we have imposed are: 1. At every time step the values of 
all relevant quantities inside a sphere of given radius cen- 
tered on the star are reset to the free flowing solution; 2. A 
uniform wind enters the computational box from the right 
hand side. We have put reflecting conditions on the axis of 
symmetry. On the left hand side of the box there is a free 
outflow boundary condition: all values of variables in the 
ghost zones are equal to the values in the corresponding 
active zones (extrapolating the flow beyond the bound- 
ary), taking care to avoid the possible inflow of material 
du e to ro ll-up features (Stone & Norman 1992 , LeVeque & 
al. 1997 ). This boundary condition, despite its simplicity, 



works quite well, and is actually used by various codes. We 
know from solar- wind simulations that the choice of such 
boundary conditions may influence the evolution in the 
case of subsonic outflow, but even in more sophisticated 
codes, the only way to solve the problem seems to take 
boundaries very far from the zone of interest (Pauls & al. 



1995). We have tested the code with our simple buondary. 



and we have found that there are no significant problems 
with refiected spurious waves, and even following the evo- 
lution of our simulations we do not find spurious wave. 
Both the stellar and the external winds are taken with 
high Mach numbers (hypersonic limit). In fact the pulsar 
wind has a thermal pressure much lower than its ram pres- 
sure and the same holds true for the ISM, because pulsars 
have typ ical p eculiar velocities of 200 kms^^(Cordes & 
Chernoff 1998 ) while the typical temperature of the ISM 
is of the order of 10* - 10^ K. 

The pulsar wind is highly re lativi stic, with a Lorentz 
facto r up to 10^ (Rees & Gunn 1973 Kennel & Coroniti 



1984 ), thus in order to handle it properly a relativistic 



code is required; however the geometry of such nebulae 
is essentially determined by the momentum flux of the 
winds (Wilkin 1996 ), so even a classical model can give a 
reasonable approximation. We decided to use a different 
adiabatic coefficient for the pulsar wind material and the 
ambient medium, so we have modified the Roe Riemann 



Solver (Roe 1981) used by the code, adding a tracer (\) 
to separate the ambient medium (with adiabatic coeffi- 
cient 7i — 5/3) from the relativistic (with 72 — 4/3) fiuid 
coming from the star (Shyue 1998, Karni 1994). We then 
integrate the set of equations: 
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where p is the density, Vz and Vr are the components of the 
velocity, P E and H are the pressure, energy density and 
enthalpy of the fluid, (p takes the value 0i for the ambient 
medium and 02 for the relativistic material. So, the value 
of the adiabatic coefficient 7 in each cell is determined by 
the relation: 
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(3) 



The physical conditions of such a bow-shock guarantee 
conservation of energy, and little mixing between the two 
media, so that we can use a fiuid treatment. However we 
added a little diffusion to avoid the growth of numerical in- 
stabilities ("carbuncle") due to the solver, and associated 
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with grid aligned shocks. In fact these instabilities develop 
on the axis and affect both the bow-shock and the contact 
discontinuity and then move towards the tail of the neb- 
ula, leading to a complete disruption of its structure. The 
introduction of a little diffusion, which avoids the growth 
of numerical instabilities, but may also damp the physical 
ones, should not be considered only as a numerical trick; 
in fact the four known nebulae show a well-defined shape, 
so there must be some process (perhaps connected to the 
magnetic field) which tends to stabilize the flow. 

2.1. Choice of initial parameters 

We have performed two different types of simulations: the 
first one limited to the head of the nebula, using a finer 
grid, and the second one extended to the tail up to about 
4 times the distance from the stagnation point to the star, 
with a coarser resolution. 

The unit length we have used, W, has the value 
lO^^cm, chosen to reproduce typical pulsar bow-shock 
nebula dimensions (Eq. The simulations of the head 
were performed on a 0.40 [/ZxO.58 W box, using a grid of 
480x700 cells with the star located in (100,0). For the 
simulations of the tail, we used a 1.45 Wx 1. 00 W box with 
a grid of 580x400 cells and the star in (420,0). The steUar 
wind has a velocity of 1000 kms~^, a density of 0.04 cm^'^ 
and a Mach number of about 13 upstream of the termina- 
tion shock along the bow shock axis. We have performed 
various simulations using different velocities for the exter- 
nal flow (400, 300, 200, 150 kms~^) and densities vary- 
ing accordingly, in order to keep the ram pressure con- 
stant. The following figures refer to the simulation of 400 
kms^^(density of 0.25 cm~^). The time needed to reach 
the steady state condition is of the order of 4-5 times the 
external flow crossing time in the computational box. 

3. Comparison with the analytic model 

3.1. Shape 

The geometry and internal structure of the bow shock 
nebula resulting from our simulations are shown in Fig. 1 
and Fig. 2, while a grey-scale contour plot of the den- 
sity is shown in Fig. 3. The positions of the termination 
shock, bow shock and contact discontinuity do not change 
among the mentioned above simulations, and also the po- 
sition of the sonic point remains essentially the same. This 
is because the solution, in the hypersonic limit, depends 
only on the ratio between ram pressures of the stellar and 
ambient medium, but not separately on their densities or 
velocities: therefore we are confident that we have approx- 
imated well the hypersonic limit. The dashed lines rep- 



resent three analytic solutions (Wilkin 1996), one scaled 
to match the contact discontinuity near the head, one to 
match the bow-shock near the head, and the long-dashed 
one, which is evaluated for the correct wind/environment 
ram pressure ratio. It can be seen that the analytic model 
reproduces quite well the shape of the nebula, although an 



uncertainty remains in the size. Fig. 1 shows that the re- 
gion of unperturbed stellar wind extends in the tail up to a 
"spherical shock" , whose position is strictly connected to 
the value of the thermal pressure in the ambient medium. 
In fact in the tail the pressure tends to the values of the 
ambient medium so the distance of the spherical shock 
can be easily derived as: 



J'back 



v//:/4^cPo = doMo7i 



1/2 



(4) 



where Mq is the external flow Mach number. 

We remark that the outflow boundaries we use may in- 
troduce spurious waves, and this plays a role in the struc- 
ture of the subsonic region (region A) especially as it may 
change the position of the triple point (where the spherical 
shock meets the termination shock). However the head of 
the simulated nebula, and the supersonic region of the ex- 
ternal shocked layer are not affected, so that the following 
comparison with the analytic model is still valid. 

Looking at the head (Fig. 2) we notice a bump in the 
contact discontinuity, due to an instability that has not 
been damped completely by diffusion. The contact dis- 
continuity in the tail (Fig. 1 and Fig. 3) tends to reach 
a limiting distance from the axis (i.e. the shocked pulsar 
wind is confined to a cylinder), with radius of the same or- 
der of the stagnation point distance; instead, the external 
bow-shock tends to reach the Mach cone. 

An upper limit to the radius of the contact discon- 
tinuity can be evaluated using an asymptotic treatment. 
Asymptotically in the tail the bow shock degenerates in 
a Mach cone, so the pressure in the external layer is es- 
sentially the same as the thermal pressure of the ambient 
medium Pq. Inside the contact discontinuity there are two 
different regions: the innermost (labeled as region A) cor- 
responds to the material shocked in the spherical shock, 
which moves subsonically; the other region (region B), 
which surrounds the previous one, is filled with the mate- 
rial which crossed the termination shock at distances from 
the pulsar less than the spherical shock. As shown in Fig. 1 
and Fig. 2, the material in the latter region moves super- 
sonically. In a steady state model the pressure in these 
two regions must asymptotically match the conditions in 
the ambient medium. So the pressure inside the contact 
discontinuity is also Pq and this justifies Eq. (^ for the 
position of the spherical shock. 

If we consider a fluid particle moving in region B, 
from a position near the stagnation point (labeled with 
subscript i„i) towards the asymptotic region in the tail 
(labeled with subscript asy), the entropy conservation re- 
quires: 



Pasyl Pjly — Pini/ Pini (^) 

On the other hand we have: Pasy = Poi and: 

Pim^ A{-fi)poV^, Pini = Pups B{'J2), (6) 

where pups is the value of the density upstream of the wind 
termination shock along the axis, and : 
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B{i)- ^m^(i)^- (7) 

If we consider an external wind with high Mach num- 
ber, the position of the spherical shock is far from the star 
so we can assume that all the material of the stellar wind 
is confined to region B and that the distance Z of the con- 
tact discontinuity from the axis is much larger than the 
radius of internal region A. The matter flux conservation 
law allows us to write: 

'^T^d^PupsVwind — VasyPasyTrZ'^ (8) 

From a two thin layer model wc found that Vasy , taken 
to be constant over the whole transverse section, is close 
to Vwind, the stellar wind velocity. This can be easily un- 
derstood because, if we assume that the tail extends to in- 
finity (this constraint will be released below), the velocity 
scale for the stellar wind material, in the absence of mix- 
ing with the ambient medium, is Koind (in the hypersonic 
limit the dependence on the external flow Mach number 
may be neglected). More specifically from the Bernoulli 
equation: 

^asy T2 ^asy '72 Pini 

2 72 - 1 pasy 72-1 pini 

together with Eq. ||, it follows: 

Vasy = V^,nd\] 1 - {71 /A(7i ) } ^ (10) 

Substituting p^sy and Vasy in Eq- (||)) using respec- 
tively Eq. (||) and Eq. ([lO|), we obtain: 

Z^^ 4d2{7iA/2A(7i)}- / 

S(72)Vl-{7iMoV^(7i)}~ (11) 

This can be considered as an upper limit to the radius 
of the contact discontinuity. For the stellar wind material 
to expand sideways beyond Z, some process is required 
that reduces the velocity in the layer or increases its pres- 
sure. This can be achieved via processes like mass accre- 
tion, or if some further shock develops. 

The distance of the bow-shock from the star given by 
the analytic model is 0.175 Ul (long-dashed line) while 
the two limiting analytic solutions of Fig. 1 and Fig. 2 
(dashed lines) are at 0.22 and 0.29 Ul. So the analytic 
model for the correct wind/environment ram pressure ra- 
tio gives a value which corresponds to the termination 
shock, while the thickness of the external layer is 0.025 
Ul, in good agreement with what was expected (Paper I). 
Even though the shocked layers of the nebula are not thin 
at all, we see that the analytic solution reproduces quite 
well the morphology, at least of the external layer, which 




2G0 500 400 500 



Fig. 1. Tail simulation. TS is the termination shock, SS is 
the spherical shock, CD is the contact discontinuity (dot- 
dashed line) , BS is the bow shock. The short-dashed lines 
represent two analytic solutions, the long-dashed line is 
the analytic solution for the correct wind/environment 
ram pressure ratio, the dotted lines indicates the sonic 
surfaces (the subsonic region is in the head). The circle 
corresponds to the region in which the stellar wind values 
are fixed every time step. 



corresponds to the Ha emitting region. In the tail simula- 
tion we see that the contact discontinuity has not reached 
the asymptotic regime, its distance from the axis being 
less than the maximum value given by Eq. (pT|). 

The behaviour of the fluid in regions A and B may 
affect the synchrotron emission from the relativistic elec- 
trons and positrons coming from the pulsar and shocked in 
the termination shock. As we can see from Fig. 1, material 
flowing in region B becomes supersonic soon in the head. 
In a relativistic case, the Lorentz factor associated with 
the bulk motion of the fluid in the tail becomes greater 
than the Lorentz factor associated with the thermal mo- 
tion of relativistic particles. If particles have relativistic 
motions their synchrotron emission can be detected only 
if their pitch angle is close to the observer direction. In 
the tail the relativistic bulk motion produces a beaming 
of the emission in the backward direction, so that the syn- 
chrotron emission from relativistic particles, in the case of 
a direction of observation perpendicular to the pulsar ve- 
locity direction, can be detected only in the head, while it 
fades in the tail. This holds true for the less energetic par- 
ticles responsible for the radio emission and may even be 
true for those responsible for the X-ray emission. In region 
A the flow is subsonic, so the radio and X-ray emission can 
be detected. In this case we expect a cylindrical emitting 
region. If the spherical shock is far from the pulsar it may 
occur that the number of high energy particles confined 
in region A is so small that their emission might be un- 
der the threshold of detection, so that this region can be 
observed only in the radio band. 
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Fig. 2. Head simulation. All curves and labels are defined 
as in Fig. 1. 



3.2. Internal magnetic field 
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Fig. 3. Tail simulation. Density contours and shades in 
logarithmic scale are shown. The little spreading of the 

contact discontinuity is due to the diffusion wc have in- 
troduced to avoid the growth of numerical instabilities. 



magnetic field wound around the symmetry axis and pas- 
sively advected by the fluid. We considered the case of 
a pulsar wind with the energy associated with the or- 
dered magnetic field much smaller than the ram pressure 
of the relativistic particles, and with the ratio between the 
two being independent of the direction. The value of such 
a field increases in the head, due to the compression in 
the termination shock, ranging from 7 times its upstream 
value, right after the shock, up to a factor of 30-35 with 
respect to its upstream value, near the contact disconti- 
nuity. In the tail it tends to decrease with respect to its 
value in the head duo to the subsequent expansion of the 
fluid until it approaches a constant value: 



7 B„. 



^Bups {7iM2^(7i)}=- 5(72)^ X 

1^^-1/4 

1 - (7iM2/^(7i)) 



(12) 



In fact the contact discontinuity tends to a limit sec- 
tion and all the magnetic field produced by the star is 
confined in such region. If the original magnetic field is 
not too far below equipartition, it can dominate the fiuid 
dynamics both in the tail and in the head. So an hydro- 
dynamical approach is asymptotically valid if: 



Pa 



» 1 



(13) 



Following a procedure analogous to the one used to 

derive the asymptotic axial distance of the contact dis- 
continuity we are able to relate the asymptotic condition 
to the initial values in the head: 



Using the numerical hydrodynamic solutions described 
above we have also evaluated the behaviour of a toroidal ^'^"v 



{7iM2A(7i)} ^= 
(7iM2/^(7i)) 



B{l2) 

1/2 



(14) 



So, even if initially the magnetic field's energy is below 

equipartition, it can reach equipartition if the external 
wind Mach number is high enough. Full MHD simulations 
are under development. 

3.3. Surface density and tangential velocity in the 

external layer 

To verify that the two thin layer model could be taken as 
a reasonable approximation of the true structure of the 
nebula, at least as far as the external layer is concerned, 
and to evaluate how far in the tail it could be used, we 
have extracted from our simulations the surface density 
and the tangential velocity in the external layer, between 
the bow shock and the contact discontinuity. 

Fig. 4 compares the surface density in the external 
layer, extracted from the tail simulation, with that evalu- 
ated by a two thin layer model. The numerical values have 
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Fig. 4. Comparison of the surface density in the external 
layer, obtained from the tail simulation (crosses), with the 
values obtained from an analytic best-fit model (continu- 
ous line) , and those from the analytic solution for the cor- 
rect wind/environment ram pressure ratio (dashed line). 



been obtained by integrating in the direction perpendic- 
ular to the bow shock. The difficulty in determining how 
good the analytic model approximation is follows from the 
fact that the analytic model depends on just one param- 
eter, while the real problem has intrinsically two length 
scales (the stagnation point distance and the thickness in 
the head). In fact we have found two limit solutions (see 
Fig. 1 and Fig. 2), so that we could put only an upper 
and a lower limit to the length scale that should be used 
in the analytic model. The analytic curve that best re- 
produces the numerical points gives an external density 
of 0.26 cm~'^ (instead of the true value 0.25 cm~^) and 
a stagnation point distance of 0.235 Ul, compatible with 
the range of the two limit solutions. The fluctuations of 
the numerical points are caused by round-off errors and 
uncertainties in the evaluation of the proper orthogonal 
direction, due to small oscillations in the shape of the bow 
shock, caused by discretisation errors in the evaluation of 
bow-shock position, eventually enhanced by the interpo- 
lation routine we use to find the orthogonal direction. In 
the simulation of the head we have noticed that, near the 
axis, the density values are lower than expected, and this 
effect corresponds to the bump of the contact discontinu- 
ity, that we have previously discussed. If corrected for the 
presence of the bump also in this zone we find reasonable 
values for the density. 

In the same way, we have evaluated the tangential ve- 
locity in the layer between the bow shock and the contact 
discontinuity. Fig. 5 shows the numerical values of the tan- 
gential velocity in the external layer, normalized to the 
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Fig. 5. Comparison of the tangential velocity in the ex- 
ternal layer, obtained from the tail simulation (crosses), 
with the velocity obtained from an analytic best-fit model 
(continuous line), and those from the analytic solution for 
the correct wind/environment ram pressure ratio (dashed 
line). 



external wind value, and compares it to the two thin layer 
model. The analytic curve that reproduces the points has 
been obtained using a value of 0.25 Ul for the distance of 
the stagnation point, a value in agreement with the two 
limit solutions. The fluctuations in the numerical points 
are much smaller than those seen in the surface density, 
which implies that the velocity pattern is more uniform in 
the layer with respect to the density. 

The fact that the analytic curves that best reproduce 
the numerical values correspond to different values of the 
stagnation point distance tells us that the density and 
velocity cannot be reproduced exactly by the same thin 
layer model. It should be stressed, though, that, in spite 
of the non-thin structure of the external layer, the thin 
layer model is not a bad approximation. Deviations from 
the continuous curves are present only in the tail where a 
steady state regime may have not been achieved yet. 

4. Penetration thickness for neutral hydrogen 

As we have demonstrated in a previous paper (Paper I) 
the penetration thickness r for neutral hydrogen in the 
external layer, depends on the ambient density and the 
pulsar velocity. For a pulsar moving in a partially neutral 
medium ( the IS M ionization fractio n is of order 0.1; Taylor 
& Cordes 1993 ; Dickey & Lockman 199C), this parameter 
plays an essential role in determining the true size of the 
nebula and the Ha photon flux. The penetration thick- 
ness depends on the process that could drag the atoms: 
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Fig. 6. Penetration thickness of the external layer for neu- 
tral hydrogen coming from the ambient medium. Uq is 
the ambient particle density, ao is the Bohr Radius. The 
four cases represent different velocity of the pulsar: A 150 
kms-i; B 200 kms'^ C 300 kms'^; D 400 kms'^ 

essentially by charge-exchange and by ionization via the 
interaction with the protons confined in the external layer. 
The cross sections of suc h proc esses are stron gly d epen- 
dent on velocity (McClure 



1966, Newman et al. 1982: Peek 



1966| ; Ptak & Stoner |1973| ), and variations of the surface 
density in the external layer give rise to changes of the 
thickness moving away from the stagnation point. In our 
previous paper we have evaluated the thickness using the 
two thin layer model limited to a region near the axis. Our 
aim is to verify if such an approximation is acceptable and 
how far in the tail can it be extended. 

The thickness of the external layer has been evaluated 
integrating the ionization and charge- exchange probabili- 
ties along the path of a neutral test hydrogen atom. So we 
have taken into account variations of density, relative ve- 
locity and temperature. We have verified that, for typical 
pulsar velocities, collisional ionization can be neglected, 
even in non equilibrium cases. 

Fig. 6 shows the penetration thickness, normalized to 
the external density. The differences in the values for the 
various simulations are due to variations of the cross sec- 
tion with velocity, while the effect of temperature is less 
sensitive. The thickness decreases at higher velocities as 
expected: in fact it follows the behaviour of the cross sec- 
tion for charge-exchange, which is the dominant process. 
Concerning the variation of t with the distance from the 
stagnation point, we can see that it has a homogeneous 
trend for different velocity values. Moving away from the 
stagnation point, the surface density increases and the 
same holds for the cross section of charge- exchange, due 



to the decrease of relative velocity between the neutral hy- 
drogen atoms coming from ISM and the protons flowing 
in the external layer. This effect, that tends to increase 
the penetration thickness, is however reduced because the 
interaction rate, which depends on the relative velocity, 
decreases even if the cross section is higher, so that the 
mean free path is longer. The thickness, at least for the 
head, ranges only over a factor 2. Thus we can say that the 
model restricted to the stagnation point, which we have 
proposed (Paper I) , can reproduce within the same factor 
the behaviour of the whole head. This can be quite im- 
portant, because, if the nebula is thin to the penetration 
of neutral hydrogen, it can remain thin in the tail, even 
far from the pulsar. 

5. Conclusion 

We have performed hydrodynamical simulations of the in- 
teraction with an ambient homogeneous medium, of the 
wind of a pulsar moving supersonically. The use of a two 
component (the ambient medium and the pulsar wind) 
code allows us to verify directly if the two thin layer model 
gives a reasonable approximation. We have also extended 
our simulations to the tail. Our results show that the an- 
alytic model is a good approximation for the shape, the 
surface density and the tangential velocity in the external 
layer, and can be extended in the tail up to about 3 or 4 
times the stagnation point distance from the pulsar. For 
the range of pulsar peculiar velocities that we have take 
into account, we have verified that the penetration thick- 
ness for neutral hydrogen in the head of the nebula has a 
rage of a factor of 2. 

This tells us that models based on the analytic two 
thin layer solution could be considered a reasonable start- 
ing point, over which we might introduce effects of mass 
loading. So it will be possible to derive a simple solution 
for the case in which the neutral component of the ambi- 
ent medium cannot be considered as a small perturbation. 
The analysis we have performed gives information on the 
asymptotic behaviour of such nebulae. As we have pointed 
out in the introduction, two out of the four known objects 
seem not to correspond to our expectations. This is evi- 
dence that in such cases we cannot assume a purely hy- 
drodynamical steady state condition, and we require some 
processes in the tail (mass loading, shocks) to model the 
observations. However the observational data on these ob- 
jects are not good enough to allow for a clear choice of the 
more adequate description. Our simulations show that the 
penetration thickness of the nebula changes slowly moving 
away from the head, thus, if we assume that ionizing pro- 
cesses other than collisions are active, mass loading may 
affect the fluid motions not only in the head but also in 
the tail. It is known that mass loading or charge-exchange 
from slowly moving particles in a supersonic fluid tends 
to reduce the velocity and to increase the pressure. The 
shocked material coming from the stellar wind becomes su- 
personic soon in the head. In the supersonic region, mass 
loading due to ionization of the neutral atoms, penetrat- 
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ing from the ambient medium, may be due to ionizing 
photons coming from the hot shocked ISM, while pho- 
toionization from the pulsar may be neglected (Paper I). 
So mass loading may justify the presence of a conical tail, 
which requires an internal thermal pressure greater than 
what we found from purely hydrodynamical simulations, 
and this suggests the necessity of ionizing processes other 
than collisions. 

We have also pointed out that the magnetic field might 
play an important role both in the head and in the tail, 
leading to a faster coUimation of the flow. The presence 
of a component of the magnetic field wounded around the 
axis of the nebula may reduce the effect of mass load- 
ing preventing the flow in region B from expanding side- 
ways and coUimating it. In fact in the case of the "Guitar 
Nebula" (the only for which we know the shape even far 
from the pulsar) we see that there is a conical tail near the 
star (which we suppose due to mass loading) which stops 
increasing and become a cylindrical tube. To justify this 
change in the shape we need some other process (that we 
believe could be the action of magnetic fields, which may 
reach equipartition) that coUimates the tail. Full MHD 
simulations are under development, as well as relativistic 
HD simulations, and will be presented in a future paper. 
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